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ABSTRACT 

We used a sample of Kepler candidate planets with orbital periods less than 200 days and radii between 
L5 and 30 Earth radii (/?©) to determine the typical dynamical spacing of neighboring planets. To derive the 
intrinsic (i.e., free of observational bias) dynamical spacing of neighboring planets, we generated populations of 
planetary systems following various dynamical spacing distributions, subjected them to synthetic observations 
by the Kepler spacecraft, and compared the properties of observed planets in our simulations with actual Kepler 
detections. We found that, on average, neighboring planets are spaced 21.1 mutual Hill radii apart with a 
standard deviation of 9.5. This dynamical spacing distribution is consistent with that of adjacent planets in 
the Solar System. To test the packed planetary systems hypothesis, the idea that all planetary systems are 
dynamically packed or filled to capacity, we determined the fraction of systems that are dynamically packed 
by performing long-term (10^ years) numerical simulations. In each simulation, we integrated a system with 
planets spaced according to our best-fit dynamical spacing distribution but containing an additional planet on 
an intermediate orbit. The fraction of simulations exhibiting signs of instability provides an approximate lower 
bound on the fraction of systems that are dynamically packed; we found that >31%, >35%, and >45% of 
2-planet, 3-planet, and 4-planet systems are dynamically packed, respectively. Such sizeable fractions suggest 
that many planetary systems are indeed filled to capacity. This feature of planetary systems is another profound 
constraint that formation and evolution models must satisfy. 



Subject headings: methods: statistical ■ 
lites: detection 



planetaiy systems - planets and satellites: general - planets and satel- 



1. INTRODUCTION 

We examine the question of whether planetary systems gen- 
erally consist of closely spaced planets in packed configu- 
rations or whether planets in the same system are generally 
more widely spaced apart. Here we adopt the traditional def- 
inition of dynamical spacing as the separation between adja- 
cent planets in terms of their mutual Hill radius, and we define 
a planetary system to be dynamically packed if the system is 
"filled to capacity", i.e., it cannot accept an additional planet 
without leading to instability. 

The degree of packing in planetary systems has impor- 
tant implications for their origin and evolution. It has 
been codified in the packed planetary systems (PPS) hy- 
pothe si s (e.g., B arnes & Quin n 2004; Raymond & Barnes 
l2005t iRavniond et alj 120061; iBames & GreenbergI 12007^ 
the idea that all planetary systems are dynamically 
packed. Previous works have invoked the PPS hypothe- 
sis to predict the existence of additional planets in sys- 
tems with observed planets locate d far apart with an in- 
termediate stabi l ity zo ne (e.g. Menou & Tabachnik 2003; 
Barnes & Oulmil l2004t "Raymond & Barnes 2005; Ji et al. 
2005t Rivera & Haghighipour 2007; Raymond etal. 2008; 
Fang & Margot^ i2012 b). since the PPS hypothesis requires 



that an undetected planet is located in that stability zone. Sys- 
tems that are observed to have dense configurations could sup- 
port the PPS hypothesis if they were shown to be dynam- 
ically packed. Such systems may i nclude KepIer-11, wit h 
six transiting planets within 0.5 AU (iLissauer et al.ll201 Ibh . 
Kepler-36, whose 2 kno wn planets have s emi-major axes 
differing by only -10% (ICarter et all 120121) . and KOI-500, 
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which has 5 planets al l within an orbital period of 10 days 
(iRagozzine etalJl2012h . 

In this study we seek to investigate the underlying distribu- 
tion of dynamical spacing in planetary systems by fitting to 
the observed properties of Kepler planet candidates. By un- 
derlying or intrinsic, we mean our best estimate of the true 
distribution of dynamical spacing between neighboring plan- 
ets in multi-planet systems, i.e., free of observational biases. 
After we derive the underlying distribution of the dynamical 
spacing between planets, we create planetary systems whose 
planets have separations that obey this distribution. We then 
subject these planetary systems to N-body integrations to ex- 
amine their stability properties, which allows us to determine 
if they are dynamically packed or not. By determining the 
fraction of systems that are packed, we can test the PPS hy- 
pothesis. 

In a related study published by Fang & Margoj (l2012al) . we 
investigated the underlying multiplicity and incUnation dis- 
tribution of planetary system s based on the Kep ler catalog 
of planetary candidates from Batalha et al.l (120121) in Febru- 
ary 2012. We created population models of planetary systems 
following different multiplicity and inclination distributions, 
simulated observations of these systems by Kepler, and com- 
pared the properties of detected planets in our simulations 
with the properties of actual Kepler planet detections. We 
used two types of observables: numbers of transiting systems 
(i.e., numbers of singly transiting systems, doubly transiting 
systems, triply transiting systems, etc.) and normalized tran- 
sit duration ratios. Within our orbital period and planet radius 
regime (P < 200 days, 1.5 < ^ < 30 /?©), we found that 
most planetary systems had 1-2 planets with typical inclina- 
tions less than 3 degrees. In the present study, we build upon 
and extend this previous investigation to explore the under- 
lying distribution of dynamical spacing in planetary systems 



using data from the Kepler mission. 

This paper is organized as follows. In Section l2n we de- 
fine our stellar and planetary parameter space. We also de- 
scribe how we created model populations of planetary sys- 
tems and how we compared them to the properties of Kepler 
planetary candidates. In Section l272l we present the best-fit 
model representing our best estimate of the intrinsic distribu- 
tion of dynamical spacing in planetary systems. In Section[3j 
we compare this distribution of dynamical spacing with that of 
the Solar System. We also make comparisons with two other 
systems, Kepler- 11 and Kepler-36, to quantify how rare such 
systems are. In Section]!] we test and quantify whether such a 
distribution of dynamical spacing implies that planetary sys- 
tems are dynamically packed, by performing an ensemble of 
N-body integrations. We briefly describe implications for the 
PPS hypothesis. Section ]5] summarizes the main conclusions 
of this study. 

2. DYNAMICAL SPACING OF PLANETS 
2.1. Methods 

Our methods for deriving the intrinsic dynamical spacing 
of planetary systems are as follows. First, we created model 
populations of planetary systems obeying different underlying 
distributions of dynamical spacing. Second, we performed 
synthetic observations of the planetary systems in these pop- 
ulations by the Kepler spacecraft. At this stage we identi- 
fied which simulated planets were detectable by the Kepler 
telescope, and which were not. Third, we compared the re- 
sulting distribution of dynamical spacing of detectable plan- 
ets from synthetic populations with that of the actual Ke- 
pler detections. The actual distribution can be easily ob- 
tained from Kepler transit data with an assumed planet radius- 
mass relationship. Mos t of these steps are fully explained in 
iFang & MargotI (I2012al) . and we refer the reader to that paper 
for details. In the following paragraphs, we summarize the 
mos t salient points of our procedure as well as any differences 
wifli lFang & Margoli (1201 2ah . 

Each model population consists of about 10^ planetary sys- 
tems, and we created various model populations that followed 
different underlying distributions of multiplicity, inclinations, 
and dynamical spacing. To generate these populations, we 
needed to restrict the range of physical and orbital proper- 
ties of the stars and planets that we considered in our sim- 
ulations. We selected ranges that would adequately over- 
lap those of a K epler sample that can be considered nearly 
complete (iHowa rd et al. 2012; Youdin 2011). Stellar proper- 
ties such as radius 7?*, stellar noise ct*, effective temperature 
Teff, surface gravity parameter \og{g), and Kepler magnitude 
K were randomly draw n from the Kepler Input Catalog (see 
iFang & Margoli]2()12a l). We only considered bright solar-like 
stars that obeyed the following ranges: 



4100 K<reff< 6100 K, 

4.0 < \og{g [cm s-2]) < 4.9, 
K<\5 mag. 



(1) 



Planet radii and orbital periods were drawn from debiased dis- 
tributions, and we obtained these debiased distributions by 
converting the observed sample of Kepler Objects of Inter- 
est (KOI; bas ed on detections up to Quarter 6 released in 
February 2012. lBatalha et al.l20I2l) into a deb iased sample us- 
ing calculations of detection efficiencies (see Fang & MargoJ 
]2012a) . We filtered the KOI sample (and correspondingly 



limited the parameter space of the synthetic populations de- 
scribed in this paper) to the following orbital period P, planet 
radius R, and signal-to-noise ratio (SNR) boundaries: 



P < 200 days, 
1.5 R(s <R<30Rq., 
SNR(-> Q8)> 11.5. 



(2) 



These limits were imposed in order to choose a sample of 
planets with properties unlikely to be missed by the Kepler 
detection pipeline. For SNR, we required an SNR>10 for 
observations up to Quarter 6, which corresponds to about 
SNR> 11.5 for observations up to Quarter 8 by assuming that 
SNR roughly scales as \/N, where is the number of ob- 
served transits. This scaling is performed because the SNRs 
of observed KOIs have be e n re ported for observations up to 
Quarter 8 in iBataIha et al.l (1201 2). whereas the actual detec- 
tions have been reported up to Quarter 6 only. Planet masses 
M were calculated by converting from planet radii R. We 
used a broken log-linear M{R) prescription ob tained by fitting 
to mas ses and radii of transiting planets (see iFang & MargotI 
]20T2al) : 
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Additionally, we repeated all of the methods described in 
this section by using an alternate mass-radius relationship: 
(M/Mgj) = (R/R(sf °^ dLissauer et all 1201 Ibl) . By adopt- 
ing this alternate mass-radius equation, our results showed 
the same best-fit dynamical spacing distribution as defined 
in Equation (]7]i with cr =14.5 (see results presented in Sec- 
tion 12.21 ). We note that both of these mass-radius equa- 
tions were obtained by fitting to the sample of planets with 
known masses and radii. Errors in the mass-radius relation- 
ships can potentially affect our results because our determi- 
nation of dynamical spacing is a direct function of plane- 
tary masses. In Figure JT] we investigate how uncertainties 
in the mass-radius relationship map into dynamical spacing 
uncertainties. Specifically, we plot histograms showing how 
the observed dynamical spacing changes if there is a l-u in- 
crease or decrease in mass for the mass-radius equation. For 
masses lower than nominal, adjacent planets appear to be less 
closely spaced and so the distribution shifts to the right. For 
masses higher than nominal, adjacent planets appear to be 
more closely spaced and so the distribution shifts to the left. 
The shifts are quantified at the end of Section l2!2] 

For orbital eccentric i ties, we adopted circular orbits, as we 
did in lFang & Margoti (l2012al) . Eccentricities do not directly 
affect our calculation of dynamical spacing, as we will define 
below in Equation (]5j. 

Regarding the multiplicity distribution in our model pop- 
ulations, we used a bounded uniform distribution with A = 
1.5-3.5 with increments of 0.25 to assign the number of plan- 
ets per system. A bounded uniform distribution has a sin- 
gle parameter A and is defined as follows: first, draw a value 
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Figure 1. (Top) Observed transiting exoplanets with known 
masses and radii shown in green, and the corresponding mass- 
radius relationship. (Bottom) Dependence of the observed 
dynamical spacing distribution on the choice of mass-radius 
relationship. In both panels, the nominal mass-radius rela- 
tionship as defined in Equations (O and (|4| is shown in black, 
and the mass-radius relationship shifted by one sigma to lower 
(higher) masses is shown in brown (red). A represents the 
number of mutual Hill radii between adjacent planets in multi- 
planet systems and is defined in Equation (|5]l. 

A^max (maximum number of planets) from a Poisson distribu- 
tion with parameter A that ignores zero values, and second, 
draw the number of plan ets from a discrete unif orm distribu- 
tion with range 1 -A^max dFang & Margotll2012al) . Thus, each 
planetary system will have at least one planet. For the inclina- 
tion distribution of planetary orbits, we used a Rayleigh dis- 
tribution with cr = 1,2° as well as a Rayleigh of Rayleigh dis- 
tribution with O'er = 1 , 2° . A Rayleigh of Rayleigh distribution 
has a single parameter and is defined as follows: first, draw 
a value cr from a Rayleigh distribution with parameter a^, and 
second, draw a value f or inchnation from a R ayleigh distribu- 
tion with parameter a (Lissau er et al.ll20"l Ibl) . These specific 
multiplicity and inclination distributions were chosen because 
they yielded fits consistent with transit num bers and transit 
duratio n ratios from Kepler detections (see iFang & Margoti 
l2012al) . Combinations of these specific multiplicity and in- 
clination distributions add up to a total of 36 possibilities. 
The difference between model populations generated in 



iFang & Margoti (l2012ah and the model populations generated 
in this study is the treatment of planetary separations, since 
here we wish to determine the underlying dynamical spacing 
of planetary systems. We used a separation criterion A to 
assess the dynamical spacing between all adjacent planets in 
multi-planet systems, where A is defined as (e.g., iGladmaiJ 
IT99llCham"bers et al.l[T996h 
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(6) 



In these equations, a is the semi-major axis, Rh 1,2 is the mu- 
tual Hill radius, and M is the mass. Subscripts *, 1, and 2 
refer to the star, the inner planet, and the outer planet, respec- 
tively. For a two-planet system not in resonance, the planets 
are required to be spaced with A > 3.46 in order to be Hill 
stable. 

In our model populations, adjacent planets in multi -planet 
systems were spaced according to a prescribed A distribution. 
We used a shifted Rayleigh distribution, which is the same as 
a regular Rayleigh distribution except shifted to the right by 
3.5 (since we require this distribution to provide values of A 
that meet the minimum Hill stability limit). Such a distribu- 
tion, as we will show, matches the observed sample well. The 
mathematical form of a shifted Rayleigh distribution / is 



/(A): 



A-3.5 



-(A-3.5)V(2<T^) 



(7) 



and is described by a single parameter a. In our model pop- 
ulations, we explored values of a =10-20 with increments of 
0.5 for a total of 21 possibilities. We chose this range of 
cr values based on the location of the observed A distribu- 
tion (blue histogram in Figure |2]i with its approximate peak 
at about 20 mutual Hill radii. This chosen range of a val- 
ues allowed us to explore different distributions of dynamical 
spacing that spanned a reasonable range of possible model A 
distributions. Increments of 0.5 were chosen as a trade-off 
between resolution and computational limitations. As will be 
seen in Section 12.21 the statistically good match between the 
data and the best-fit model demonstrates that our increments 
are sufficiently small and have appropriately sampled the pos- 
sible range of A distributions. 

In order to enforce that adjacent planets are spaced ac- 
cording to the prescribed A distribution, we performed the 
following steps. For each synthetic planetary system, the 
first planet's orbital period is drawn from the debiased pe- 
riod distribution. If the system's multiplicity is greater than 
one, for the second planet we draw its separation A from the 
first planet using the prescribed A distribution and we also 
draw a value from the debiased period distribution. If that 
value is less/greater than the first planet's period, then the sec- 
ond planet will be the inner/outer planet and its exact period 
will be calculated by satisfying the drawn A separation from 
the first planet. This process repeats if the system has addi- 
tional planets. These steps are different from Fang & Margo^ 
(|2012a), where in that study all planetary orbital periods were 
chosen by drawing them from a debiased distribution. We 
verified that the periods drawn to match the A distribution 
provide a very close match to the debiased period distribution 
as well. 



3 



After the creation of each model population, we performed 
synthetic observations of each population's planetary systems 
by the Kepler spacecraft in order to determine which planets 
were transiting and detectable (see Fang & Margot 2012a). 
The transiting requirement was evaluated by picking a ran- 
dom line-of-sight (i.e., picking a random point on the celestial 
sphere) and computing the planet-star distance projected on 
the plane of the sky. The minimum of that distance was com- 
pared to the radius of the host star to determine if the planet 
in our simulations transited or not. The detection requirement 
was assessed by calculating each transiting planet's SNR, de- 
fined as 



(8) 



where the first fraction gives the depth of the transit, repre- 
sents the number of transits up to Quarter 6, and cr, represents 
stellar noise (Combined Differential Photometric Precision 
or CDPP; Christiansen et al. 2012). Since CDPP is quarter- 
to-quarter dependent, we used the median CDPP value over 
all available quarters. In the calculation of SNR, we took 
into account gaps between Kepler quarters, the fact that not 
all sta rs are observed each quarter, and a 95% duty cycle 
jFang & Margot 2012a). If the calculated SNR for a transit- 
ing planet met or exceeded the SNR threshold for detection 
(SNR=10), then it was considered detectable. 

Lastly, we determined the goodness-of-fit between each 
model population's detected planets and the actual Kepler de- 
tections. This was ascertained by comparing the A distribu- 
tions of adjacent planets in their multi-planet systems. We 
performed a Kolmogorov-Smirnov (K-S) test to assess the fit 
between the A distributions, and this comparison yielded a 
p-value that we used to evaluate the null hypothesis that the 
distributions emanate from the same parent distribution. This 
K-S probability was used to determine how well a particu- 
lar model matched the observations. We also calculated the 
goodness-of-fit for multiplicity (by comparing with observed 
Kepler numbers of transiting systems using a chi-square test) 
and for inclination (by comparing with observed, normalized 
transit duration ratios using a K-S test) to check that th ey were 
consistent with the data (see iFang & Margolll2bl2al) . While 
we only generated model populations with underlying multi- 
plicity and inclination distributions that are considered to be 
good fits to the data based on our previous work, this extra 
step allowed us to confirm that any models with acceptable A 
fits also produced acceptable multiplicity and inclination fits 
to the Kepler data. We determined which model populations 
were most consistent with the data by combining (multiply- 
ing) the probabilities associated with each one of the 3 statis- 
tical tests that probed multiplicity, inclination, and dynamical 
spacing. We assumed that these probabilities are independent. 

Accounting for all combinations of multiplicity, inclination, 
and dynamical spacing distributions, in total we created 756 
model populations with about 10^ planets each. As described 
earlier, each of these model populations underwent synthetic 
observations by Kepler as well as statistical tests. The next 
section presents our results. 

2.2. Results 

We report our results on the dynamical spacing (represented 
by the criterion A) in multi-planet systems based on Kepler 
data. Using the methods described in the previous section, we 
determine that our best-fit model for the intrinsic A distribu- 



tion is a shifted Rayleigh distribution (see Equation d?)) with 
0-= 14.5. 

This best-fit distribution is plotted in Figure |2] where we 
show its probability density distribution (top panel) as well as 
its cumulative probability distribution (bottom panel). This 
best-fit distribution has a mean value of A = 21.7 with a 
standard deviation of 9.5. About 50% of neighboring planet 
pairs have A separations larger than 20, and about 90% of 
neighboring planet pairs have A separations larger than 10. 
This best-fit distribution was obtained by considering shifted 
Rayleigh distributions with increments in a of 0.5. The 
mean values for distributions with a = 14.0 and a = 15.0 are 
A = 22.3 and A = 21 .0, respectively. The combined probabil- 
ities for the match to the data using distributions with a = 14.0 
and a= 15.0 are less than half the combined probability using 
0-= 14.5. 

Our results are valid for the range of stellar and planetary 
parameters given in Equations ([T) and (|2), most notably a 
minimum planet radius of 1 .5 and a maximum orbital pe- 
riod of 200 days. It is possible that planets are actually even 
more closely spaced than this best-fit distribution if there are 
planets located in intermediate locations with radii less than 
1.5 /?0. Therefore, our findings about the A distribution can 
be used to represent the spacing of planetary systems confined 
to the scope of our study, or can serve as an upper limit for the 
spacing of planetary systems that include planets with smaller 
radii. 

Figure [3] shows the comparison between the A distribution 
of simulated detections from this best-fit model's population 
and the observed A from actual Kepler detections. Note that 
this figure does not show the underlying A distribution that 
is plotted in Figure |2j instead, this figure shows the distribu- 
tion of simulated planets that would have been detected, in 
order to make an appropriate comparison with the observed 
A distribution. The K-S test for comparing these two dis- 
tributions yields a p-value of 56%, indicating that we cannot 
reject the null hypothesis that these distributions are drawn 
from the same parent distribution. In other words, this model 
is consistent with the observations. 

Comparison between Figures|2]and|3]shows that the under- 
lying A distribution is similar to the observed A distribution- 
both distributions have peaks near A ^ 20. This suggests that 
the observed A distribution is not indicative of a significant 
population of non-transiting and/or low-SNR planets (within 
the planet radius and orbital period limits of our study) lo- 
cated in-between detected planets; otherwise, the underlying 
A distribution would have on average lower A values than 
those of the observed distribution. We caution again that our 
study cannot rule out the existence of a population of planets 
with R<\ .5R(^, so the actual underlying A distribution could 
be different from what our results indicate. The similarity be- 
tween the observed and underlying A distributions is due to 
the rarity of cases where a system has Aobsei-ved > Aunderiying 
(i.e., an undetected planet located in-between two detected 
planets). The stringent geometric probability of transit means 
that outer planets are more easily missed than inner planets. 
As a result, it will be rare to find cases where an intermedi- 
ate planet is non-transiting and therefore missed, but both the 
innermost and outermost planets are transiting and detected. 
For such a case, which rarely occurs, the observed A would 
be greater than the underlying A. 

We discuss how we expect these results to change if an al- 
ternate mass-radius relationship is used. In particular, Fig- 
ure[T]shows how the observed dynamical spacing distribution 
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Figure 2. The best-fit model's underlying A distribution is shown as a magenta-colored curve. This distribution represents our best 
estimate of the true or intrinsic (i.e., free of observational bias) distribution of dynamical spacing between all neighboring planets 
meeting our orbital period and planet radius criteria. Recall that A represents the difference in semi-major axes between two 
adjacent orbits; it is expressed in units of the mutual Hill radius. The best-fit model is a shifted Rayleigh distribution as defined 
in Equation dTji with a = 14.5; the top plot depicts the probability density and the bottom plot shows the cumulative probability. 
The black-colored curves show the range and sampling frequency of model distributions that follow different cr parameter values 
ranging from <t = 10 to cr = 20 with increments of 0.5, as examined in this study. 



changes depending on various choices for the mass-radius re- 
lationship. For computational expediency, we chose to evalu- 
ate the effects of the mass-radius relationship on the observed 
A distribution. We use this as an approximation for the effects 
on the underlying A distribution, with the justification that the 
observed and underlying A distributions appear similar (see 
previous paragraph). For the mass-radius relationship shifted 
down by 1 sigma, the best-fitting shifted Rayleigh distribution 
has (7 = 16.5. For the mass-radius relationship shifted up by 
1 sigma, the best-fitting shifted Rayleigh distribution has a 
= 12.5. As a result, we estimate that our derived dynamical 
spacing distribution can span a =12.5-16.5 due to uncertain- 
ties in the mass-radius scaling. 



3. COMPARISON TO THE SOLAR SYSTEM, KEPLER-1 1, AND 
KEPLER-36 



We can compare our results (i.e., the intrinsic dynamical 
spacing or A distribution) with the dynamical spacing distri- 
bution of the Solar System, if we extrapolate beyond the ra- 
dius and period limits (Equation ^) of our study. Figure |4] 
shows our intrinsic A distribution of planetary systems over- 
plotted with a histogram of the A distribution of the Solar 
System planets. From this figure, it is interesting to note that 
the distributions appear to be relatively similar and that the 
Solar System planets may be similarly spaced as most exo- 
planets in general. A K-S test between our cumulative A dis- 
tribution and the sample of A values between adjacent planets 
in the Solar System yields a p-value of 66.2%, indicating that 
the Solar System A distribution is consistent with that of Fig- 
ure |2] 

The orbital evolution of the planets i n the Solar System is 
known to be chaotic and unstable (e.g.,|S ussman & Wisdom 
[T988t lLaskari[T989[ [T990t iSussman & Wi sdom 1992; Laskm 
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Figure 3. The top plot shows the comparison between the ob- 
served A distribution from actual Kepler detections (shown in 
blue) and the A distribution from simulated detections in our 
best-fit model population (shown in orange); the K-S proba- 
bility for this match is 56%. The bottom plot shows this com- 
parison in cumulative form and its histogram points are con- 
nected by lines for easy viewing, and also shows the compar- 
ison with the properties of detected planets from other model 
populations (a =10 to a =20 with increments of 0.5; shown 
in black). The 99.5% confidence region of acceptable fits in- 
cludes model populations with a ranging from 12.5 to 17.0. 



119941; iMichtchenko & Ferraz-Menoll2001l; iLecar et al.ll200n) . 
The inner Solar System can be potentially unstable within 
the Sun's remaining lifetime due to a secular reso- 
nance (Batvsin & Lauahlin 2008). Laska] (11994 and 
iLaskar & Gastineau (2009.) have shown that inner planets 
could be ejected or collide. Numerical simulations of the 
planets in the outer Solar System suggest that they are 
packed (Barnes & Quinn 2004; Raymond & Barnes 2005; 
iBarnes et al.ll2008l) . All of these results suggest that the So- 
lar System is dynamically packed. If we consider the Solar 
System to be dynamically packed then it is possible that other 
planetary systems with similar planet multiplicities and A dis- 
tributions are also dynamically packed. This prompted us to 
verify whether planetary systems in general are dynamically 
packed (see Section|4]l. 

Kepler- 11 is a planetary system with six known transit- 
ing pla nets in a closely spaced configuration (iLissauer et al.l 
1201 lal) . All six transiting planets have orbits smaller than the 
orbit of Venus, and five of the six transiting planets have or- 
bits smaller than the orbit of Mercury. This appears to be 



Figure 4. Comparison of the cumulative A distribution be- 
tween the best-fit model's A distribution (i.e., shifted 
Rayleigh distribution with a = 14.5; purple solid line) and the 
Solar System's A distribution (i.e., histogram based on its 8 
planets; dotted green line). 

a very compact system, and we calculate the A separations 
of the innermost five planets to be Ab-c =5.7, Ac-d =14.6, 
Ad-e =8.0, and Ag-f =11.2. We did not calculate the A sep- 
aration of the f-g planet pair because the mass of planet 
g is not known and only has an upper limit . Accounting 
for the 1(7 uncertainties on mass reported bv ILissauer et al.l 
(1201 lal) . the dynamical spacing of these pairs have the follow- 
ing ranges: Ab-c =5.1-7.0, Ac-d =13.0-17.3, Ad-e =7.2-8.8, 
and Ag-f =10.0-12.6. All of the Kepler-1 1 planets are within 
the planet radius and orbital period scope of our study, and we 
apply our knowledge of the intrinsic dynamical spacing (i.e., 
Figure|2|i to this system. We find that the separation Ab-c =5.7 
is more closely spaced than 98.9% of adjacent planet pairs 
in multi-planet systems, Ac-d =14.6 is more closely spaced 
than 74.6%, Ad-e =8.0 is more closely spaced than 95.3%, 
and Ae-f =11.2 is more closely spaced than 86.8%. These 
high percentages indicate that the planetary separations in the 
Kepler- 11 system are much smaller than average separations 
in planetary systems, and we conclude that Kepler- 11 is un- 
usual in terms of the density of its configuration. 

Kepler-36 has two known transiting planets with a large 
density contrast (their densities differ by a factor of ^8) yet 
they orb it closely to one a nother (semi-major axes differ by 
--10%) dCarter et all 1201 2 j). Such close orbits with dissimi- 
lar densities are unusual compared to the planets in the Solar 
System, where the denser terrestrial planets are located in the 
inner region and the less-dense giant planets are located in the 
outer region. We calculate the dynamical spacing between 
the two planets in Kepler-36 to be A=4.7. In comparison to 
our intrinsic A distribution of dynamical spacing, a separation 
of A=4.7 is more closely spaced than 99.7% of neighboring 
planet pairs of planetary systems in general. 

4. DYNAMICAL PACKEDNESS OF PLANETS 

Section lZ2l described the best-fit A distribution of planetary 
systems based on Kepler data, with a mean value of A = 21 .7 
(Figure |2]l. In this section, we investigate whether this distri- 
bution of A implies that planetary systems are dynamically 
packed or not. By dynamically packed, we refer to a plan- 
etary system that is filled to capacity and cannot include an 
additional planet without leading to instability. 

To investigate whether planetary systems are dynamically 
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packed, we performed long-term N-body integrations of plan- 
etary sys tems gen erated for our best-fit model population (see 
Sections 12. 1112.21 1. For each multiplicity (i.e., 2-planet sys- 
tems, 3-planet systems, 4-planet systems), we randomly chose 
1000 planetary systems for which we performed long-term in- 
tegrations. In total, we performed 3000 integrations. We did 
not include single planet systems because they are irrelevant 
for studies of dynamical packedness and we did not include 
systems with multiplicities higher than 4 plan ets because they 
are rel atively rare for our parameter space (iFang & MargotI 
l2()T2ah . 

For each integration, we added an additional planet when 
testing for stability; this planet had a mass equal to the low- 
est mass of all original planets and its initial conditions in- 
cluded an orbital eccentricity of zero, an inclination drawn 
from a Ray leigh of Rayleigh distribution with <Jc = 1° (see 
Section ITl] ). and random values for its argument of pericen- 
ter, longitude of the ascending node, and mean anomaly. This 
additional planet was placed in-between the orbits of exist- 
ing planets, and if there were 3 or 4 original planets, we ran- 
domly determined which 2 adjacent planets would be receiv- 
ing a new neighbor The additional planet's semi-major axis 
was calculated so that it was located with equal mutual Hill 
radii distances between its neighboring planets. These initial 
conditions for the additional planet (e.g., low eccentricities, 
low inclinations, a mass equal to the lowest mass of origi- 
nal planets, a semi-major axis located at equal A distances 
from neighboring planets) are very conservative in the sense 
that we have chosen initial conditions that are very amenable 
to stability, as we determine whether a planetary system with 
this additional planet can remain stable or not. 

Our simulations were performed using a hybrid 
symplecti c/Buhrsch-Stoer integrator from the Mercury 
package (I ChambersI Il999h . and we used a timestep that 
covered 1/25 of the innermost planet's orbital period. Simula- 
tions were performed for a length of 10^ years; the instability 
timescales had median values less than 10^ years. Possible 
outcomes included either a stable system with no instabilities 
or a system with at least one instability defined as a collision 
between the star and a planet, a collision between planets, 
and/or an ejected planet. All of these planetary systems were 
verified to be stable for 10** years before adding the additional 
planet. 

The results of our simulations can be divided into two 
camps. The first group is composed of planetary systems that 
became unstable in our integrations. This suggests that these 
planetary systems are dynamically packed, since the addition 
of another planet in an intermediate orbit resulted in an unsta- 
ble planetary system. We found that 31%, 35%, and 45% of 
2-, 3-, and 4-planet systems were unstable, respectively (Ta- 
ble [U. The second group is composed of planetary systems 
that did not exhibit any signs of instability. For these systems, 
although they were stable within the scope of our integrations, 
they may still be dynamically packed. Possible reasons in- 
clude: an instability may occur on a longer timescale than our 
integration time, there may be additional planets in the sys- 
tem beyond the scope of our orbital period range of 200 days, 
or there may be additional planets in the system smaller than 
1.5 7?0, which is the minimum radius of our study. All of 
these factors would affect the determination of the stability of 
the system, and so for this second group of systems we are 
agnostic about their dynamical packedness. 

Accordingly, we can only confidently provide a lower limit 



Table 1 

Lower Limits on the Percentage of Dynamically Packed Systems 



System Multiplicity 


Percentage of Packed Systems 


2-Planet Systems 


> 31% 


3-Planet Systems 


> 35% 


4-Planet Systems 


> 45% 



Lower limits on the percentage of dynamically packed systems as 
obtained from the fraction of numeiical integrations exhibiting instabil- 
ities. A planetary system is considered to be dynamically packed if the 
addition of another planet causes instability. The results of our simula- 
tions only provide lower limits because the absence of instability does 
not indicate that a system is not dynamically packed (see main text). 

on packed systems by concluding that at least 31-45% (de- 
pending on the system's multiplicity) of planetary systems 
with dynamical spacings consistent with our best-fit A distri- 
bution (Figure|2]i are dynamically packed. These lower limits 
are also presented in Table [T] Note that systerns with lower 
multiplicity are more common (IFang & Margotll2012al) . 

The packed planetary systems (PPS) hypothesis is the idea 
that all planetary systems are dynamically packed, and there- 
fore cannot hold additional planets without becoming unsta- 
ble. The results of our long-term numerical integrations are 
consistent with the PPS hypothesis, as we find a sizeable 
lower limit of 31-45% (depending on the system's multiplic- 
ity) of planetary systems to be dynamically packed. 

5. CONCLUSIONS 

We have generated model populations of planetary systems 
and simulated observations of them by the Kepler spacecraft. 
By comparing the properties of detected planets in our simula- 
tions with the actual Kepler planet detections, we have deter- 
mined the best-fit distribution of dynamical spacing between 
neighboring planets. This best-fit distribution is our best esti- 
mate of the underlying (i.e., free of observational bias) dis- 
tribution of dynamical spacing for our orbital period {P < 
200 days) and planet radius (1 .5 < R <3Q R^) parameter 
regime. Stemming from this distribution, the main results of 
this study are: 

1. On average, neighboring planets are spaced 21.7 mutual 
Hill radii apart, with a standard deviation of 9.5. This distance 
represents the typical dynamical spacing of neighboring plan- 
ets for all systems included in the parameter space described 
above. 

2. Our best-fit distribution of dynamical spacing is consis- 
tent with the dynamical spacing of neighboring planets in the 
Solar System, with a K-S p-value of 66.2%. If we consider the 
Solar System to be dynamically packed, then it is not unrea- 
sonable to ask whether other planetary systems with similar 
dynamical spacing are also dynamically packed. 

3. Based on our best-fit distribution of planetary spacing: 
>31% of 2-planet systems, >35% of 3-planet systems, and 
>45% of 4-planet systems are dynamically packed. This 
means that such systems are filled to capacity and cannot hold 
another planet in an intermediate orbit without becoming un- 
stable. 

4. Our results on the dynamical packedness of planetary 
systems are consistent with the packed planetary systems hy- 
pothesis that all planetary systems are filled to capacity, as we 
find sizeable lower limits on the fraction of systems that are 
dynamically packed. 

5. Compact systems such as Kepler-11 and Kepler-36 rep- 
resent extremes in the dynamical spacing distribution. For ex- 
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ample, the two known planets in Kepler-36 are more closely 
spaced than 99.7% of all neighboring planets represented by 
our orbital period and planet radius regime. 

We thank the referee for valuable comments. This research 
was performed using resources provided by the Open Sci- 
ence Grid (OSG), which is supported by the National Science 
Foundation and the U.S. Department of Energy's Office of 
Science. 
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